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We propose a flat-histogram Monte Carlo method to efficiently sample fractal landscapes such as 
escape time functions of open chaotic systems. This is achieved by using a random- walk step which 
depends on the height of the landscape via the largest Lyapunov exponent of the associated chaotic 
system. By generalizing the Wang-Landau algorithm, we obtain a method which simultaneously 
constructs the density of states (escape time distribution) and the correct step-length distribution. 
As a result, averages are obtained in polynomial computational time, a dramatic improvement 
over the exponential scaling of traditional uniform sampling. Our results are not limited by the 
dimensionality of the phase space and are confirmed numerically for dimensions as large as 30. 
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The development of Monte Carlo methods had a dra- 
matic impact on our understanding of high-dimensional 
systems. The spectrum of applications of these meth- 
ods was considerably expanded with the development of 
optimized methods, such as flat-histogram [US] and par- 
allel tempering [4 , and now-a-days includes problems in 
a variety of fields, ranging from fluid dynamics [5j and 
spin systems jTHS] to protein simulations O [7]. These 
methods efficiently compute averages using non-uniform 
sampling and are optimized to problems on which the 
high dimensionality of the system leads to phase spaces 
with rough and complex energy landscapes. 

In chaotic dynamical systems, complex landscapes ap- 
pear even in low dimensions due to the sensitivity of ini- 
tial conditions. Prominent examples of such landscapes 
appear in open systems showing chaotic transients. Tran- 
sient chaos is a classical problem of nonlinear dynam- 
ics [5] with recent applications in fields ranging from 
quantum scattering to chemical and biological reactions 
in fluid flows [HUTU]. In open chaotic systems, the number 
of trajectories with escape time t decays as p(t) ~ e~ Kt (k 
is the escape rate) and the t = oo trajectories build a un- 
countable fractal set. The dependence of t on the phase- 
space variables r build thus a fractal landscape as the 
one illustrated in Fig. [T] This example of extreme rough 
landscape poses major numerical challenges [9]- While al- 
gorithms beyond uniform sampling have been proposed 
for specific problems, e.g. to compute the fractal dimen- 
sion [TT] or to find long-living trajectories [P2rtl4| . there 
is still no general framework to sample the phase space 
of such systems. 

In this paper we show how Monte Carlo methods can 
be applied to fractal landscapes such as those appearing 
in open chaotic systems. The crucial step is to design 
a random walk able to explore the extreme roughness 
of fractal landscapes. We show that an efficient flat- 
histogram simulation is only obtained using a random- 




FIG. 1: (Color online) Fractal landscapes appear in the 
escape time function of open chaotic systems. Escape time 
t as a function of the phase space coordinates (2/1,2/2) at 
xi — X2 = of the 4-Dimensional coupled Henon maps de- 
fined in Eq. jjj). Inset: magnification of the red rectangle 
showing a Cantor-set-like profile. The circles (states) and ar- 
rows (proposals) represent the random walk (green/red indi- 
cate accepted/rejected proposals) underlying the Monte Carlo 
sampling. 

walk step length a which scales with the landscape height 
t as a(t) ~ e~ Xht , where Xr, is the maximum Lya- 
punov exponent of the underlying chaotic system. More- 
over, by extending the Wang-Landau procedure [5] to 
the proposal distribution, we obtain an adaptive algo- 
rithm which provides simultaneously p(t) and <r(t). In 
open chaotic systems, our approach changes the scaling 
of the computational effort from exponential to polyno- 
mial with the maximum escape time and both efficiently 
finds the large t trajectories and computes averages over 
the phase space. 

We consider a fractal landscape as an escape time func- 
tion of a transient chaotic system. Given a discrete-time 
open dynamical system r t+1 = F(r t ) defined in a D di- 
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mensional phase space f2, the escape time t(r) is defined 
as the number of iterations needed for an initial condi- 
tion r to leave the region of non-trivial dynamics [9]. 
We propose an algorithm that constructs both the total 
volume p(t) of the landscape and the optimal proposal 
width a(t) at each i, in a predetermined time-spectrum 
[iminjimax] an d with a precision / which is successively 
reduced (initially f — e and a(t) = p(t) = 1 for all t). We 
consider a state dependent random walk in the space of 
initial conditions r G [25 initialized at r £ T,t = t(r), 
which evolves according to the following four steps: 

51- propose a state r'eT with t' — i(r') G [imim t max ] 
[e.g., using Eq. ([I]) below]. 

52- accept/reject the state according to flat-histogram 
choice [Eq. ^ below]. 

53- update p(t) and a(t) as: 

53.1- p(t) <- p(t)f (Wang-Landau); 

53. 2- cr(t) <- a(t)f if t 1 = t; <r(t) <- a{t)/f if t' < t. 

54- After a number of repetitions of S1-S3, refine / -s— 
V7 and go to SI. 

After the convergence (/ = / m i n > 1), the random walk 
corresponds to a flat-histogram Monte Carlo simulation 
on t Q] (using only SI and S2). We now describe in more 
detail the steps S1-S4. 

Si-Proposal - The ideal proposal should be able to ex- 
plore the order of the landscape for an efficient search. In 
discrete spaces, often considered in spin systems, there 
is a natural local proposal given by flipping a single 
spin [TS]. In continuous spaces the locality of the pro- 
posal is determined by the neighborhood around the 
present state. Fractal landscapes do not have a global 
characteristic length scale [SI [9] and therefore we consider 
a time dependent proposal width a — <x(i). Accordingly, 
we choose the conditional probability of proposing a new 
state r' given r as 
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where u is a random vector on the surface of the D- 
dimensional unit-hypersphere, and a(t) gives the charac- 
teristic length of the distribution [26] . 

We now show how a(t) has to scale with t for an ef- 
ficient proposal. We consider the construction of the 
Cantor set [SI H] as a paradigm of fractal landscape, see 
Fig. [2] The construction starts by splitting the interval 
[0, 1] in the intervals [0, 1/a], [1/a, 1-1/6], [1 - 1/6, 1] 
and assigning the escape time t = to the middle in- 
terval (plateau at t — 0). This procedure is repeated 
on each of the 2 surviving intervals by assigning t = 1 
to each of their 2 middle intervals (plateaus at t = 1), 
and again in the remaining intervals ad infinitum. In 




FIG. 2: (Color online) The construction of the 2-scale Can- 
tor set with scales 1/a and 1/6 as a paradigmatic fractal land- 
scape. A plateau at t with width e(t) generates two plateaus 
at t + 1, with widths e(t + l) = e{t)/a and e(t + l) = e(t) /b, as 
indicated in the inset. At each t, the 2* plateaus have lengths 
e k (t) with k = 0...t. 



order to achieve an efficient proposal we have to know 
the scaling of the typical length of the plateaus e(t) 
with t. For the 1-scale Cantor set (a — 6), each of 
the 2* plateaus have a unique length given by e(t) = 
(l-2/a)(l/a)* and thus e(t) = e(t). For the 2-scale Can- 
tor set (a / 6), the 2* plateaus have t + 1 different lengths 
e k (t) = (l-l/a-l/6)(l/a)*- fc (l/6) fc with k = 0, ...,t and 
the number of plateaus with size e k (t) is the binomial co- 
efficient B{t, k), see inset of Fig. [2] The total length at 
t is p(t) = (1 - 1/a - l/b)(l/a + 1/6)* - exp(-«tf). The 
conditional probability of being at a plateau of length 
Ek(t) at a given t is 



P(k\t) 



P(k,t) B(t,k)e k (t) 



P{t) 
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The characteristic plateau size is thus naturally chosen as 
e(t) = Ek'(t) where k = k* maximizes P{k\t) in Eq. 
Using Stirling's approximation we obtain k* « t/(l+b/a) 
and thus 



= £k*(t) = cxp(-< 



alog(6) + 61og(a) 



(3) 



In the context transient chaos, the construction of the 
Cantor set corresponds exactly to the escape time func- 
tion of the one-dimensional open tent map [27] and the 
exponent Xl — al ° s ^^ los t a ) corresponds to its positive 
Lyapunov exponent [9 a . This leads to a natural interpre- 
tation for a choice of a(t) with e(t) as given in Eq. (jjjj: 
in order to ensure that two chaotic trajectories (initiated 
at r and r') remain correlated up to time t, their initial 
distance \r — r'\ should be reduced exponentially with 
t, with an exponent equal to the positive Lyapunov ex- 
ponent responsible for the divergence in forward time. 
In a generic fractal landscape, generated by a higher di- 
mensional system, this divergence is dominated by the 
maximal Lyapunov exponent Al and therefore 

a(t) ~ e(t) ~ e-^ (4) 

should be used in any isotropic proposal such as Eq. ([I]). 
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S2 - Acceptance - Because of the extreme roughness of 
fractal landscapes, we use a flat-histogram simulation [1 
on the variable t, which plays the role traditionally played 
by energy. The detailed balance of this Monte Carlo 
process is fulfilled when the acceptance, the conditional 
probability of accepting a proposed state r' given r, fol- 
lows the Metropolis' choice |15 



A{r 



min < 1 



p{t{r')) g{r' -> r) 
M p(t(r)) g(r -> r') 



(5) 



where g(r — > r') is given by Eqs. ([I]) and Q. Since we 
are considering projections in i, it is useful to define the 
conditional probability A{t) of accepting a proposal given 
a time t |15j . In the spirit of flat-histogram simulations, 
a signature of an efficient random walk is a A(t) which 
does not strongly depends on t. In Fig. [3] we show that 
only when the scaling in Eq. Q is used in the Eq. ([lj, 
we obtain a constant A{t) on a flat-histogram simulation 
in the generic 2-scale Cantor set and thus an efficient 
simulation. 




FIG. 3: Characteristic random-walk step a has to scale as the 
typical plateau size e in order to achieve a constant acceptance 
ratio in time. Acceptance ratio A(t) of a flat-histogram simu- 
lation on the 2-scale Cantor set with parameters (a, b) = (3, 4) 
and three (see legend) slightly different exponents A on the 
proposal width a(t) ~ e~ xt , with Xl given in Eq. For 
A < Xl and f > 1, a(t) ^> i(t), and A(t) decays (exponen- 
tially) as a consequence of lack of proposals to t' > t. For 
A > Xl and i > 1, a(t) <C i(t), A(t) increases to 1 but the 
simulation gets stuck in the same plateau as all proposals are 
for t' = t. 



S3 - Wang-Landau update - In systems on which p{t) 
and a{t) (or are known, we use steps S1-S2 to sam- 
ple them. However, for generic landscapes, p(t) and a(t) 
are unknown. We take advantage of the analogy between 
p(t) and a density of states and apply the Wang-Landau 
procedure to compute it [2]. This is done by successive 
approximating p(t) in steps S3 and S4 of our approach. 
To compute alt), we propose the following generaliza- 
tion of the Wang-Landau procedure (step S3.1) to the 
proposal distribution (step S3. 2): if the proposed state 
has an escape time smaller than the present state, t' < t, 
we decrease a(t) by dividing it by /. If it has the same 
escape time, t' = t, we increase a(t) by multiplying it by 
/. This converges asymptotically (/ — > 1) to a Markov 
process. 



S4 - Refinement - In order to efficiently reduce the 
precision parameter / we take its square-root (as in 
Wang-Landau algorithm [2]) after a number of round- 
trips [151 I16| . defined as the movement in the time- 
spectrum from t min to i max and back to t mm , and we 
use an equivalent procedure to the one in Ref. [17] in 
order to avoid error saturation. 

We now confirm the generality of the approach de- 
scribed above through numerical simulations in generic 
fractal landscapes generated by a family of N coupled 
Henon maps on a ring 



Ai - xl + By l +k(xi- Xj+i) 



(6) 



for i = 1, N, N + 1 = 1, parameters k = 0.4, B = 0.3, 
A x = 3 (if N > 1), A N = 5, and Ai = A\ + {A N - 
Ai)(i — 1)/ {N — 1). This choice of parameters ensures 
that a chaotic map is obtained in the N = 1 case and 
the map considered in Ref. [13 is recovered for N = 2 
(used as a representative case to illustrate our algorithm) . 
Initial conditions start on a 2N hypercube T = [—4, 4] 2N 
and escape is defined as leaving T. In Fig. [4] we confirm 
the convergence and validity of our algorithm by showing 
that p(t) coincides with uniform sampling, a(t) scales 
with the Lyapunov exponent reported in Ref. [13] . and 
both the acceptance and the histogram of visits to escape 
time t are flat in t. We observed similar results for all 
N's up to N = 5. 




FIG. 4: (Color Online) Confirmation that our method yields 
the correct values of k and Az, and converges to a flat his- 
togram simulation for the case N = 2 in Eq. (6}. (a) p(t) 
obtained through our method (line) and through uniform 
sampling (circles) with the exact same computational effort. 
Lower inset: histogram H(t) of visits to escape times t. (b) 
a(t) obtained through our method. The dashed line shows the 
scaling c _Al ' with Xl ~ 1-33 obtained in Ref. [13] . Lower in- 
set: the acceptance ratio A(t). We used [t m m,imax] = [1,45], 
log e / m in = 2 -13 and all plotted quantities were measured on 
the last refinement. 

We now compare our approach to uniform sampling 
in terms of computational efficiency. For each i max , we 
compute the average number of map iterations n(i max ) 
per sampled state with t = i max . This comparison guar- 
antees that the maximum uncertainty of any observable 
to be measured on the landscape is the same in both 
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cases for all t < i max . For a uniform-sampling simula- 
tion, n(t max ) ~ l/p(i m ax) ~ e Ktmax . For a flat-histogram 
simulation, obtained after the convergence of our method 
S1-S4, we adopt a conservative approach which avoids 
the sampling of correlated states by considering a single 
sample of i max for each round-trip. The estimation of 
n{t max ) in this case is based on the number of steps per 
round-trip expected of an unbiased random walk with lo- 
cal steps (At « 1), which scales as ~ i^ax- Additionally, 
each proposal requires t map iterations and, since the 
histogram is flat, for each round trip one gets an addi- 
tional £ max contribution, leading to an expected scaling 
of n ~ ^max- Figure [5] confirms the dramatic improve- 
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FIG. 5: The computation effort in number of map itera- 
tions n(t max ) scales polynomially with maximum escape time, 
in contrast to the exponential scaling of uniform sampling. 
Main panel: results for uniform sampling (solid line) and flat- 
histogram simulation achieved by our method (squares) for 
the representative case N = 4 in Eq. (6|. The dashed line 
indicates the scaling t^ax- Inset: dependence of the efficiency 
on the phase-space dimension 2N at a fixed i max = 16. 

ment from exponential (uniform sampling) to polynomial 
(our approach) scaling in the coupled Hcnon maps. The 
significance of these results become apparent by notic- 
ing that i max = 237 (last point in Fig. [5]) corresponds 
to p(t max ) rj 1CJ~ 109 , meaning that we are able to sam- 
ple extremely rare states. For such level of accuracy, our 
method requires an implementation with arbitrary pre- 
cision [T8] which in our case was able to resolve states 
which differ by KT 137 (since cr(i max = 237) « KT 137 ). 
Interestingly, the slight but clear deviation from the pre- 
diction t^ax seen in Fig. [5] shows that flat-histogram sim- 
ulations on fractal landscapes are not purely diffusive on 
t, a phenomenon known in spin-systems as critical slow- 
ing down [TTU [20] ■ This phenomenon is enhanced with 
increasing dimension and contributes to the exponential 
increase of n max with N for a fixed i max , as shown in the 
inset of Fig. [5j Still, an uniform sampling in such a high 
phase-space dimension (27V = 30) would need impracti- 
cable n ~ 10 34 map iterations to sample one state with 

^max — 16. 



In summary, we have shown how flat-histogram Monte 
Carlo simulations can be performed on fractal land- 
scapes. The crucial ingredient is to consider a random- 
walk step size dependent on the height of the landscape. 
The correct dependency should scale as the characteristic 
length of the landscape and can be obtained through an 
adaptive procedure which generalizes Wang-Landau's al- 
gorithm to the proposal distribution. This idea can find 
applications in any rough landscape for which there is a 
typical width at a fixed height and a strong dependence 
with the height. Fractality can be considered as an ex- 
treme case of roughness which naturally occurs in the 
escape-time landscape of open chaotic systems. In this 
case, our results show that the Lyapunov exponent Xl, 
a fundamental property of the chaotic dynamics, is an 
essential ingredient for a flat-histogram simulation. 

We emphasize the significance of our results for numer- 
ical investigations of open chaotic systems. Our method 
automatically provides k and Al, is not limited to low 
dimension, and allows for the computation of expected 
values of any observable using a flat-histogram simula- 
tion. For the specific problem of finding the chaotic sad- 
dle |12H14j . which is indirectly solved in our simulations 
by storing trajectories with large t, our findings show that 
best results are achieved using a proposal which scales as 

More generally, besides high-dimensionality, the sensi- 
tivity of initial conditions in chaotic systems is a major 
reason for using statistical methods in physics. Monte 
Carlos methods in dynamical systems were traditionally 
limited to uniform sampling and, only recently, opti- 
mized methods (with non-uniform sampling) were ap- 
plied for the problem of finding trajectories with low 
chaoticity [31] [52]. Our approach opens the perspec- 
tive of using the full strength of optimized Monte Carlo 
methods to problems which involve the computation of 
averages in chaotic systems. Spatially extended [53] and 
non-hyperbolic Hamiltonian 24J systems are natural can- 
didates for future applications of this approach. 
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